template<class Resample,class Kernel>
class Ghkc
{
	public:
	Ghkc(Resample *R, Kernel *K){

		boost::random::uniform_01<> U;
		_runif=new RandomG::Random<boost::random::uniform_01<> >(U);
		_R=R;
		_K=K;
		d=new boost::math::normal_distribution<>(0,1); 
		rt=new rtnormNC();
		rt2=new rtnorm();
		boost::random::normal_distribution<> N2;
		_rnorm= new RandomG::Random<boost::random::normal_distribution<> >(N2);
		sigma=0.000000001;
	}
	~Ghkc()
	{
		delete d;
		delete _runif;
	}

	void ghk_oneside(double *L, double* trunpt, int *above, int *dim,double *eta, int *n, double *res)
	/*	routine to implement ghk with a region defined by truncation only on one- side
							r mcculloch 8/04
		if above=1, then we truncate component i from above at point trunpt[i-1]
		L is lower triangular root of Sigma
		random vector is assumed to have zero mean
		n is number of draws to use in GHK	
		modified 6/05 by rossi to check arg into qnorm
		modified 2012 by ridgway to make it awesome!!!! 
	 */
	{
	   int i,j,k;
	   double mu,tpz,u,pa,pb,arg;
	   mat z(*n,*dim);
	   double *prod=new double[*n];
	   Normalize(prod,n);
	   double thres=*eta*(*n);
	   z.fill(0);
	   double *Nt=new double[(*n)*(*dim)];
	  // sobol::Owenscramble O;
	   //O(*n,*dim,0,Nt);
	   *res = 1;
	   int kk=0;
	   int b=0;
	 // mat X=vect2mat<double>(L,*dim,*dim);
	 //cout << X;
	 //mat U(1,*dim);
	 //U.fill(0);
	//SampleGibbs(dim,above,trunpt,L,&U); 
	      for(j=0;j<*dim;j++) {
			//Resample
		      double ess=Ess(prod,n);
		      if(ess<thres){	
		      cout << "ess " << ess <<"\n";
				*res*=mean2(prod,n);
				//cout << mean2(prod,n); 
		      		fResample(prod,n,&z);
				b=1;
			//	_K->Move2(&(z),j);
				//cout << "m1 "<< mean(z,0);		
				for(int ll=0;ll<3;ll++)
				{
					GibbsMove2(dim,n,above,trunpt,L,&z,j);
					//RandomWalk(dim,n,above,trunpt,L,&z,j);
				}
				//cout << "m2 "<< mean(z,0);		
			}
			// set weigths	      
				for(i=0;i<*n;i++) {
					 mu=0.0;for(k=0;k<j;k++) mu += L[j*(*dim)+k]*z(i,k);
					 tpz = (trunpt[j]-mu)/L[j*(*dim)+j];
					 if(above[j]){
					    pa=0.0; pb = cdf(*d,tpz);   
					 }else{
					    pb=1.0; pa = cdf(*d,tpz);
					 }
					 prod[i]*= pb-pa;
					 //cout << prod[i] << " ";
					 u = (*_runif)(); 
					 //u=Nt[kk];
					 arg=u*pb+(1.-u)*pa;
					 if(arg > .999999999) arg=.999999999;
					 if(arg < .0000000001) arg=.0000000001;
					 z(i,j) = quantile(*d,arg);
			      		kk++;
				}
		//		cout << "\n\\\\\\\\\\\\\\\\\\\\\\\n";
	   }
		      *res*=mean2(prod,n);
	   delete[] prod;
	}

	inline void Normalize(double *w,int *n)
	{
		for(int i=0;i<*n;i++)
		{
			w[i]=1;
		}
	}
	inline double mean2(double *w, int *n){
		double sum=0;
		for(int i=0;i<*n;i++){
			sum+=w[i];
		}
		return sum/(*n);
	}
		
	inline void fResample(double *w, int *n,mat *z)
	{
		mat v=growingvect(*n);
		(*_R)(n,w,v);
		Arangemat(*z,v);
		Normalize(w,n);
	}
	inline double Ess(double *w, int *n)
	{
		double sum=0;
		for(int i=0;i<*n;i++)
		{
			sum+=w[i];
		}
		double sum2=0;
		for(int i=0;i<*n;i++)
		{
			sum2+=pow(w[i]/sum,2);
		}
		return (double)1/sum2;

	}
	inline void GibbsMove2(int *dim,int *n,int *above,double *trunpt,double *L, mat *z, int j)
	{
		int k;
		double mu,tpz,u,pa,pb,arg;
		double ltrunc=-numeric_limits<double>::infinity();
		double rtrunc=+numeric_limits<double>::infinity();
		double curr=0;
		mat Z(*n,*dim);
		Z=*z;
		int b1=1;
		int b2=1;
		for(int i=0;i<*n;i++)
		{
			for(int l=0;l<j;l++)
			{
				b1=1;
				b2=1;
				//calcul du max
					ltrunc=-numeric_limits<double>::infinity();
					rtrunc=numeric_limits<double>::infinity();
					for(int ii=l;ii<j;ii++)
					{

							double temp=trunpt[ii];
						//	cout << "trunc" << trunpt[ii] << " ";
							for(int jj=0;jj<=ii;jj++)
							{
								if(jj!=l)
								{
									temp-=L[ii*(*dim)+jj]*(Z)(i,jj);
						//			cout << (*z)(i,jj) << " ";
								}
							}
						//	cout << "\n";
							temp=temp/L[ii*(*dim)+l];
							//cout <<  "temp "<<temp << " ";
							if(L[ii*(*dim)+l]<0)
							{
							//	cout << "r" ;
								if(rtrunc>temp)
								{
									b1=0;
									rtrunc=temp;
									//cout << "rtrunc "<<  temp << "\n";
								}
							}else{
							//	cout << "l" ;
								if(ltrunc<temp)
								{
									b2=0;
									ltrunc=temp;
									//cout << "ltrunc "<<  temp << "\n";
								}
							}

					}
			//	cout << "l "<<ltrunc;
			//	cout << " r "<< rtrunc<< "\n";
			//	double me=as_scalar(sum((*z)(i,span::all))-(*z)(i,l));
				if(b1)
				{
					Z(i,l)=(*rt)(0,ltrunc,pINF,1);
				}else if(b2){
					Z(i,l)=(*rt)(0,-pINF,rtrunc,1);
				}else{
					Z(i,l)=(*rt2)(0,ltrunc,rtrunc,1);
				}
				//cout << "z " << Z(i,l);
			//	cout << "\n";
			}
			//				cout << "j "<<j<< " ";
							//cout << *z;
		}
		*z=Z;
	}
	

	inline mat Proposal(mat eta, double sigma, double* L,int t,int *dim)
	{
		mat xi(t,1);
		double sum=0;
		for(int i=0;i<t;i++)
		{
			xi(i,0)=(*_rnorm)();
		}
	/*	for(int i=0;i<t;i++)
		{
			sum=0;
			for(int j=0;j<=i;j++)
			{
				sum+=L[j*(*dim)+i]*xi(j,0);
			}
			xi(i,0)=sum;
		}
		*/
		return sigma*xi.t()+eta;
	}
	inline mat Thres(double *L,mat Y, int t, int *dim){
		mat xi(t,1);
		double sum;
		for(int i=0;i<t;i++)
		{
			sum=0;
			for(int j=0;j<=i;j++)
			{
				sum+=L[i*(*dim)+j]*Y(j,0);
			}
			xi(i,0)=sum;
		}
		return xi;
	}
	inline void RandomWalk(int *dim,int *n,int *above,double *trunpt,double *L, mat *z, int j)
	{
		int k;
		double mu,tpz,u,pa,pb,arg;
		double mu2,tpz2,u2,pa2,pb2;
		double curr=0;
		mat Y(*n,j);
		sigma/=2;
		bool t=false;
		double ratio=0;
		for(int i=0;i<*n;i++)
		{
			Y(i,span::all)=Proposal((*z)(i,span(0,j-1)),sigma,L,j,dim);	
			t=false;
			double prod=0;
			double prod2=0;
			//check constraints
			mat thres=Thres(L,Y,j,dim);
			for(int k=0;k<j;k++)
			{
				if(thres(k,0)<trunpt[k])
				{
					t=true;				
				}				
			}
			//calculate acceptance ratio
			if(!t)
			{
	      			for(int jj=0;jj<j;jj++) {
					mu=0.0;
					mu2=0;
					for(int k=0;k<jj;k++) 
					{
						mu += L[jj*(*dim)+k]*(*z)(i,k);
						mu2 += L[jj*(*dim)+k]*Y(i,k);
					}
					tpz = (trunpt[jj]-mu)/L[jj*(*dim)+jj];
					tpz2 = (trunpt[jj]-mu2)/L[jj*(*dim)+jj];
					if(above[jj]){
					    pa=0.0; pb = cdf(*d,tpz);   
					    pa2=0.0; pb2 = cdf(*d,tpz2);   
					}else{
					    pb=1.0; pa = cdf(*d,tpz);
					    pb2=1.0; pa2 = cdf(*d,tpz2);
					}
					if(pa==1)
					{
						pa=0.999999;
					}
					if(pa2==1)
					{
						pa=0.999999;
					}
					//cout << "ratio "<< pdf(*d,Y(i,jj));
					prod+=log( (pb-pa)*pdf(*d,Y(i,jj)));
					prod2+=log( (pb2-pa2)*pdf(*d,(*z)(i,jj)));
				}
				if(log((*_runif)())<(prod-prod2))
				{
					ratio+=1;
					(*z)(i,span(0,j-1))=Y(i,span::all);
				}
			//		ratio+=1;
			}
	
		}
		cout <<"ratio "<< ratio <<"\n";


	}

private:
	RandomG::Random<boost::random::uniform_01<> > *_runif;
	rtnormNC *rt;
	rtnorm *rt2;
	Resample *_R;
	boost::math::normal_distribution<> *d;
	RandomG::Random<boost::random::normal_distribution<> > *_rnorm;
	Kernel *_K;
	double sigma;
};
